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I. INTRODUCTION 

Experimental processes such as deep-inelastic scat- 
tering (DIS), semi-inclusive deep-inelastic scatter- 
ing (SIDIS) and Drcll-Yan (DY) have provided invalu- 
able information about the structure of the nucleon 
[TTj . With several new experimental facilities with 100% 
duty factor under construction, SIDIS will play an in- 
creasingly important role in the development of our the- 
oretical and experimental understanding of the structure 
of the nucleon. The elusive s — s asymmetry |12H15j 
is one area of interest that may finally be pinned down 
through the results obtained at these new facilities. The 
distribution of the spin of the proton |16H36j is an area 
of current excitement where polarized SIDIS is poten- 
tially extremely valuable through the study of tranverse 
momentum dependent parton distribution functions |37l - 
147] . which will complement work on generalized parton 
distributions [12S3HS3]. 

To allow these studies to fulfil their potential, we 
must develop a deep understanding of the fragmentation 
functions [54] . particularly their flavor, spin and trans- 
verse momentum dependence. Fragmentation functions 
appear in certain scattering reactions, for example, in 
SIDIS experiments |55l 156) and in e + e~ annihilation re- 
actions [STlI61| . Experiments are planned to use SIDIS 
to probe the flavor dependence of the parton distribu- 
tion functions (PDFs), for example, and therefore under- 
standing fragmentation functions has become very im- 
portant. Phenomenological extraction of fragmentation 
functions suffers from significant uncertainty, even for fa- 
vored fragmentation functions, which effects the system- 
atic errors associated with extracting the flavor depen- 
dence of PDFs through SIDIS. The increasing interest in 
SIDIS experiments led to the development of the NJL- 
jet model [62f[65j . which builds on the Field- Feynman 
quark-jet model (FFQJM) [55] , by using an effective chi- 
ral quark model to provide a unified framework in which 
calculations of both quark distribution and fragmenta- 
tion functions can be performed. NJL-jet model cal- 
culations of pion fragmentation functions were obtained 



in Ref. [52]. The NJL-jet model was extended to in- 
clude strange quark contributions and kaon fragmenta- 
tion functions were calculated in Ref. 63] • Further exten- 
sions of the model involved the inclusion of vector meson, 
nucleon and antinucleon fragmentation channels |64) . as 
well as the study of their transverse momentum depen- 
dence [65] and Collins fragmentation functions [67U69] . 

The probability of a fragmenting quark to produce 
two hadrons is represented by dihadron fragmentation 
functions (DFFs). DFFs have been studied recently in 
Refs. [70] [7T] in order to understand their dependence on 
invariant mass of the two produced hadrons. The focus 
of Ref. [70] was to fit parameters for a spectator model 
to output from the PYTHIA event generator [75] tuned 
for HERMES experiments [73] for DFFs with a depen- 
dence on the sum of the light-cone momentum fractions 
of the two produced hadrons and their invariant mass 
squared. Ref. [7T] focused on studying DFFs for large in- 
variant mass. DFFs with no invariant mass dependence 
were studied in the NJL-jet model in Ref. [74] at the 
model momentum scale of Qq = 0.2 GeV 2 . In order to 
compare the results with experimental data, we need to 
evolve the DFFs up to a typical experimental scale. The 
evolution equations for the DFFs are derived in Ref. [58] 
from factorization of the cross-section for the production 
of two hadrons in e + e~ annihilation in the MS factor- 
ization scheme. In Ref. [75], the non- quark evolution 
equations for DFFs were studied, while Ref. [75] focused 
on the QCD evolution equations for singlet quark and 
gluon DFFs. The ratio of the dihadron and single hadron 
fragmentation functions, which is useful when consider- 
ing experimental measurements, was also examined in 
Refs. [75] [75] . Initial conditions for DFFs for different 
pairs of hadrons and different values of z\ and Z2 are in- 
vestigated in Ref. [77], with a focus on the correlation 
function R cor obtained in the FFQJM [55] , 

An area of current interest in which the dihadron 
fragmentation functions of quarks may be useful are 
transversity distributions |31j . Transversity distributions 
are one of the three leading-twist distribution functions 
that don't vanish when integrated over the transverse 
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momentum. They describe the quark structure of the 
nucleon (the other two being unpolarized and helicity 
quark distribution functions) and these functions enter 
into asymmetries with chiral-odd versions of a special 
type of DFF known as interference fragmentation func- 
tions (IFFs) [78H82] . IFFs are DFFs with a depen- 
dence on the polarization of the fragmenting quark. In 
Refs. 83-85^, it was suggested that DFFs may be use- 
ful in extracting transversity distributions by considering 
the SIDIS production of two hadrons with small invari- 
ant mass. Transversity distribution functions are not a 
focus of this paper, but are presented as motivation for 
further investigation into DFFs. 



This work focuses on performing QCD evolution of 
the DFFs from the NJL-jet model momentum scale of 
Qq = 0.2 GeV 2 to a typical experimental momentum 
scale of Q 2 — A GeV 2 . In Section [n] we present a 
brief summary of fragmentation function equations 
from which the model scale solutions were obtained 
and used as input for the evolution equations of the 
DFFs. Section |III| describes the method for solving 
the evolution equations for single hadron fragmentation 
functions (SFFs), which are needed for the evolution 
of the DFFs. It also serves as a simple version of the 
method used to solve the DFF evolution equations, 
while the method for solving the evolution equations 
for the DFFs is described in Section |IV| A comparison 
of the model scale and evolved scale DFFs is presented 
in Section [V] Section I VII shows how the evolution code 
works on data from Ref. |76j as well as comparing our 
solutions to that data. Our data is evolved to a range of 
values of Q 2 in this section to display how the up quark 
and gluon DFFs change for larger values of Q 2 . 



II. SINGLE HADRON AND DIHADRON 
FRAGMENTATION FUNCTIONS FROM THE 
NJL-JET MODEL 

In Ref. [73], integral equations for the single hadron 
and dihadron fragmentation functions from the NJL-jet 
model are described, and the method employed to solve 
them at the model scale of Qq = 0.2 GeV 2 is presented. 
SFFs appear in the cross section for SIDIS experiments 
and thus play an important part in the theoretical under- 
standing of these experiments. In the NJL-jet model the 
SFFs, Dq(z), which correspond to the probability of pro- 
ducing a hadron h with light-cone momentum fraction z 
from a fragmenting quark q, are given by [6 2) 

The first term on the right hand side of Eq. ([!]) is the 
renormalized elementary quark fragmentation function, 
which corresponds to the process where the detected 
hadron is the only emitted hadron. We refer to this 
term as the driving function. The second term corre- 
sponds to the probability of emitting a hadron after the 
first emission step in the quark cascade and these terms 
have a sizeable effect at low values of z, while vanish- 
ing for higher z values. To solve the second term we use 
d®{z) = dg(l — z)\ h=q Q to write all functions in terms of 
their relation to the emitted hadron h. 

Dihadron fragmentation functions are another impor- 
tant tool in the theoretical understanding of the structure 
of hadrons. In the NJL-jet model, the DFF are given by 
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where the first term corresponds to the probability of 
producing hadron hi from the quark q at the first emis- 
sion step in the cascade, followed by hadron h 2 produced 
either directly afterwards or further down in the quark 
decay chain, while the second term is similar to the first 
one, except for hi -H- h 2 ■ These two terms constitute the 
driving function of the DFFs, similar to the first term 
in Eq. (JlJ. The third term on the right hand side of 
Eq. ([2| corresponds to the probability of having both the 
detected hadrons produced after the first hadron emis- 



sion. DFFs correspond to the probability of producing 
two hadrons, hi and h 2 , in the decay chain of a fragment- 
ing quark q, with light-cone momentum fractions Zi and 
Z2, respectively. 

Results for the SFFs and DFFs from the NJL-jet model 
at the model scale of Qq = 0.2 GeV 2 are described in 
detail in Ref. [74]. In this paper, they are used as the 
input for the DFF evolution equations that will be dis- 



cussed in Sections III and IV In Fig. 1(a) we present a 
3-dimensional plot of * (zi,z 2 ), at the model scale, 
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FIG. 1. 7r + 7r dihadron fragmentation function fo r th e u 

quark at the (a) model scale (Qq = 0.2 GeV 2 ) 
evolved scale (cjp = 4 GeV 2 ) 



and 



(b) 



the 



while in Fig. 1 (b) the result for the same DFF evolved to 
4 GeV 2 is shown. These plots demonstrate the effect of 
evolution on the DFFs, particularly where the functions 
achieve their peaks with respect to Z\ and z%. 



whereas the non-singlet quark fragmentation functions 
are decoupled and can be solved separately. The non- 
singlet (d^_(z,Q 2 )) and plus-type (D h q+ (z,Q 2 )^j quark 

fragmentation functions are, respectively, constructed 
from the combinations of SFFs 



D*^(z,Q 2 ) = D'l(z,Q 2 ) - D\{zM 2 ) 
= D h Jz,Q 2 )-Dl(z,Q 2 ), 



(3) 



and 



D h qt (z,Q 2 )=D h qt (z,Q 



= D^(z,Q 2 ) + Dl(z,Q 2 ), 



(4) 

where qi is the fragmenting quark. These combinations, 
rewritten using charge symmetry, allow for a simpler 
method of solving the evolution equations. 
We define the variable t as 



t=- T \n 
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where 



a s {Q 2 ) =4tt/ A, In 
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is the leading-order strong coupling constant, /3 = (33 — 
2n/)/3 is the one-loop /? function, rif is the number of 
flavors and Aqcd is the QCD scale parameteiQ We 
write the evolution equations with respect to t rather 
than InQ 2 to simplify the numerical calculation. 

The QCD evolution equations for the SFFs allow us to 
determine the SFFs at momentum scales that vary from 
the scale at which they are originally defined. This is 
achieved by calculating the rate of change of the SFF with 
respect to the momentum scale. The non-singlet, plus- 
type and gluon leading-order (LO) evolution equations 
are, respectively, given by 



^-D h -(z,t) = 
dt i, v ' > 



E 
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(7) 



III. EVOLUTION OF THE SFFS 

To evolve the DFFs, we need to evolve the SFFs as 
well. This section will focus on the evolution of the SFFs, 
and will also serve as a simple introduction to the method 
used to solve the DFF evolution equations. This proce- 
dure for solving the SFF and DFF evolution equations 
can, of course, be used for models other than the NJL- 
jet model. 

The single hadron fragmentation function evolution 
equations used in our calculations were based on those 
presented in Ref. [SS]- The evolution equations are writ- 
ten in the form of non-singlet quark, plus-type quark 
and gluon fragmentation function equations. The plus- 
type quark and gluon fragmentation functions are cou- 
pled and therefore need to be solved simultaneously, 



^-D h + (z,t) = I 
dt it { ' ' J z 



1 dy 
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In this work we take rif = 3 and Aqqjj = 0.25. 
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The left hand sides of Eqs. ([7])-(|9| represent the rate of 
change of the corresponding SFFs with respect to t. The 
right hand sides of these equations represent the effect 
that a parton j (either a quark of flavor qj or a gluon 
<?), that emits a hadron h with light-cone momentum 
fraction z/y, has on the evolution of the non-singlet (q~), 
plus-type (qf) or gluon (g) SFFs, through the splitting 
functions Pji(y) (obtained from Ref. [86 ), where i is the 
parton for the corresponding SFF on the left hand side. 



To solve Eqs. ([7|-(|9]), we express the derivatives as fi- 



nite differences using 



df(t) _ f(t j+1 )-f(tj) 



dt 



At 



(10) 



where f(t) is the corresponding SFF. We divide the range 
of interest for t into N t steps of size At. 

The integrals on the right hand side of the LO evo- 
lution equations are converted into sums over logarith- 
mically disretized values of y (denoted by zi). The 
corresponding equations for the non-singlet, plus- type 
and gluon fragmentation functions are, respectively, re- 
arranged to obtain the functions at the (k + l)th step in 
t such that 
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D h q . (z m ,t k+1 ) = D h q .(z m) t k )+At^2J2^ i P v « fr) Kj I f ■ ' 

j l=m 
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D h +{z m ,t k+1 ) = D h +[z m ,t k ) + AtJ2 
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D h Jz m ,t k+1 ) = D h (z m ,t k ) + At 



^A Zl 
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The first term on the right sides of Eqs. 



(11|)-(|13|) are 

The 



the fragmentation functions at the (fc)th step in t 
second term on the right hand side of each equation is the 
change in the fragmentation function from the (/c)th step 
to the (fc + l)th step in t. The SFF at Q\ are inserted 
as the input at k = 1, with the evolution to the next 
step, t 2 = ti + At, calculated using the previous result. 
This process is repeated to obtain the SFF evolved to the 
chosen Q 2 at fjv t +i.. 



IV. EVOLUTION OF THE DFFS 

The DFF evolution equations are derived from fac- 
torization of the cross-section for the production of two 
hadrons in e + e~ annihilation in the MS factorization 
scheme in Ref. [58]. Using jet-calculus, Ref. [87] deduces 
the evolution equations for DFFs with an explicit depen- 
dence on the invariant mass of the hadron pairs, Mh, 
which are addressed as extended dihadron fragmentation 
functions. The latter are important as they relate to 
experimental results that include the dependence on in- 
variant mass spectra. We concentrate on the DFF that 
have been integrated over the invariant mass. The LO 
evolution equation for DFFs, from Ref. [87], reads 



D^( Zl ,z 2 ,Q*) 
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where Q 2 is the momentum scale, a s (Q 2 ) is the strong and a sum over the repeated indices is implied. The rate 
coupling constant at the corresponding momentum scale 
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at which the DFFs change with respect to InQ 2 is rep- 



resented on the left hand side of Eq. (14). The first term 



on the right hand side of the LO DFF evolution equation 
represents the effect that a parton j fragmenting into two 
hadrons, hi and h 2 , has on the fragmentation of parton 
i into the two hadrons, through the splitting function 
Pji(u). The second term represents the effect of parton 
i splitting into two partons, j and k, that fragment sep- 



arately to produce hi and h 2 with light-cone momentum 
fractions u and 1 — u, respectively, through the splitting 
function P^ .(u). For the QCD evolution of the DFFs, 

both Pji (u) and P£ • (it) were obtained from Ref . [57] . 



In Eq. (14), the parton i can be either a quark, an- 



tiquark or gluon. We choose to express the evolution 
equations for the quark and gluon DFFs, respectively, 
written in terms of t (Eq. (|5|) as 



-D h q ^{z u z 2 ,t) 



du 



21+22 



hih 2 ( 
1j \ 



1-Z2 



u(l-u) 



Z-2 



1-u 1 



du 



d_ 



D h ^(z 1 ,z 2 ,t) 



1 



zi+z 2 



u 
du 
v? 



hih 2 ( 
1] \ 



1-22 



21 

1-22 



du 



-D 



hi ( z _l 



u(l-u) 93 \u' J 93 \l 



du 



u z J \u u J J z u[l — u) y 

I 



To obtain non-singlet (^D^ h2 (zi, z 2 ,t)j and plus- 
type (^D h ^ h2 (zi, z 2 ,t)j quark DFFs we use the combi- 



nations 



D h ^ (z u z 2 ,t) = D h q l h * (z u z 2 , t) - D h q l h2 (zi,z 2l t) 

= D^ h *(zi,z 2 ,t) - D^{zi,z 2l t), (17) 



and 



D 1 ^' 12 (z U Z 2 ,t) = D h q l h2 (Zi ,Z 2 ,t) + D'^ h2 ( Zl ,Z2,t) 

= D% h *{z 1 ,Z2,t) + D%f h *(zi,Z2,t), (18) 

respectively. The combination of terms on the second 
line of each equation has been rewritten using charge 
symmetry and this is the form that is employed to solve 
the LO DFF evolution equations. 

Using Eqs. (17) and (18), we write the evolution equa- 



tions in terms of the non-singlet quark, plus-type quark 
and gluon DFFs as 
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For clarity, we show the sums over the repeated indices 
and use Eq. (21 1 to display how the combinations are 



is dominated by the first term of Eq. Q._Thc model 



applied to simplify the equations. The non-singlet quark 
evolution equation is decoupled from the plus-type quark 
and gluon DFFs and can be evolved separately from 
them. Using Eq. ( 10 1 and converting integrals into sums 
over logarithmically discretized values of u, expressions 
for the DFFs evolved to the (k + l)th step in t can be 



obtained, producing results analogous to Eqs. ( 11 )-( 13 ) 



V. RESULTS 

In this section we present the results comparing the 
model scale DFFs with those evolved to Q 2 — 4 GeV 2 
for u — > 7r + ir~, u — > n + K~ and q K + K~ , where 
q = u,d, s. The first subsection explores the evolution 
of by comparing the model and evolved DFFs at 

particular values of Z\ or z%, while the second subsec- 
tion focuses on favored and unfavored hadron emission 
in the evolution of D 7 ^ K . Finally, the last subsection 



demonstrates the evolution of Dt 



for q = u, d or s. 



A. Q 2 evolution of Dl * 

We consider the DFF for an up quark fragmenting 
to 7r + and 7r~. When the up quark fragments to 7r + , 
for which it is the favored emission channel, it produces 
a down quark, which has the favored emission channel 
to 7r~. Since both emissions are favored channels for 
the detected hadrons in this quark cascade, the DFF 
has sizeable peaks in the higher z 2 an d z\ regions for 
zi = 0.5 (Fig. |2(a)[ ) and z 2 = 0.5 (Fig. |2(b)[ ), respectively. 



For Z)J 71 , the second term of Eq. Q is zero (because 
d,™ = 0) and the integral term is small, so this DFF 



scale plot for n fixed at Z\ — 0.5 (Fig. [2(a) I has the 
shape of a favored single hadron fragmentation function 
since fixing z\ effectively makes the first term on the right 
hand side of Eq. |2]) a constant multiplied by t he fav ored 
fragmentation . For z 2 fixed at 0.5 (Fig 2(b) I, the 



model scale D^ w is shaped by the elementary quark 
fragmentation function d* + , resulting in a peak at higher 
z±, while having a very small contribution at low values 
of zi. After evolution of the DFF, there is a reduction 
in magnitude and a shift in the peak towards the low 
z 2 region for zi = 0.5 (Fig. 2(a)). When z 2 is fixed at 
0.5 (Fig. 2(b)[ ), the magnitude of the DFF is reduced and 
the peak value shifts towards the low zi region. Both 
plots in Figs. [2] display a range of values at low z where 
the evolved DFF obtains a larger magnitude than the 
model scale DFF. At higher momentum scales, the low 
Zi and z 2 regions of the DFFs grow in magnitude because 
they can access the gluon emission channel. 

We present the results for z\ and z 2 fixed to 0.2 in 



Figs. 3(a) and |3(b) respectively, to investigate the DFF 
at low fixed light-cone momentum fraction. The struc- 
ture ofthe model scale D^ T for z\ = 0.2, shown in 
Fig. 3(a) is similar in shape to that of the model scale 
DFF at zi — 0.5, having the peak in the higher z 2 , ex- 
cept it is spread out more and has a lower peak value. At 
z 2 = 0.2, the structure of the model scale DFF is again 
similar to the corresponding z 2 = 0.5 plot in Fig. |2(b)| 
being small in the low zi region and having a large peak 
in the higher zi region, which is rather narrow. Evolu- 
tion of the DFF results in a shift of the peak towards 
the lower z regions, with the magnitude of the evolved 
D* 77 becoming larger than the model scale D^ T at 
mid-range values of the allowed light-cone momentum 
fraction; rather than in the lower range of values that 
was observed for the Zi and z 2 fixed to 0.5 results. The 
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FIG. 2. 7r + 7r _ dihadron fragmentation functions for the u 
quark at the model scale (Qo = 0.2 GeV 2 , shown by dotted 
red line) and the evolved scale (Q 2 — 4 GeV 2 , shown by solid 
black line) for (a) z\ = 0.5 and (b) 22 = 0.5. 



FIG. 3. tt + tv~ dihadron fragmentation functions for the u 
quark at the model scale (Qq — 0.2 GeV 2 , shown by dotted 
red line) and the evolved scale (Q 2 — 4 GeV 2 , shown by solid 
black line) for (a) zi = 0.2 and (b) 22 = 0.2. 



shape of the evolved n for 22 = 0.2 in Fig. 3(b) 
appears very similar to that of the evolved D™ +7! for 
22 = 0.5 in Fig. 2(b) [ whereas the shape for the evolved 
£>£ +7r ~ for z\ = 0.2 (Fig. 3(b) I is quite diffe rent to the 
corresponding result at zi — 0.5 in Fig. 2(a) Instead of 
the concave structure at z\ = 0.5 shown in Fig. 2(a) at 
z\ = 0.2 (Fig. |3(a) ) the evolved DFF has a large contri- 
bution at low Z2 and steadily decreases as 2 2 increases. 



B. Q 2 evolution of Dl +K 

In Figs. [4] we present the results for D* +K , where 
the up quark is a favored channel for 7r + emission, but 
the remnant down quark is an unfavored channel for K~ 

k shows no contri- 



emission. At the model scale, 
bution in th e larg e z 2 and Z\ regions for z\ (Fig. 4(a) ) 



and 2 2 (Fig. 4(b) ) fixed at 0.5, respectively. For K 
at the model scale the second term of Eq. ([2| is zero (be- 



cause d 1 ^ = 0) and the integral term is small, so 
is dominated by the first term of Eq. ([2]). In Fig. 4(a) 



the model scale DFF has the structure of the unfa- 
vored , while also being suppressed in magnitude 

by <fj (z\ = 0.5), which achieves its peak value in the 
high z\ region while vanishing in the low z\ region. For 
2 2 = 0.5 (Fig |4(b)) , the model scale DFF shows a very 
small magnitude for values of Z\ because of the combi- 
nation of multiplied by . Elementary fragmen- 
tation functions for favored emission channels are very 
small in the low 2 region, and achieve large peak val- 
ues in the high 2 region. This forces D^ +K (21,22) to 
have a very small magnitude in the low z\ region as it 
is dependent on d 7 ^ (21). is an unfavored SFF and 

therefore is constructed by the integral term on the right 
hand side of Eq. ([I]) because the first term equals zero. 
Unfavored SFFs peak in the low 2 region and have very 
small magnitudes in the medium to high 2 region. Both 
of these effects combine to cause the resultant low peak 
in the middle of the allowed region of Z\. 
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= 4 GeV 2 





(a) 2l =0.5: z 2 D* 
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(b)z 2 =0.5: z 1 Dl +K 



0.5 



Evolved Scale 
= 4 GeV 2 




0.5 



FIG. 4. ir + K~ dihadron fragmentation functions for the u 
quark at the model scale (Qq = 0.2 GeV 2 , shown by dotted 
red line) and the evolved scale (Q 2 — 4 GeV 2 , shown by solid 
black line) for (a) zi = 0.5 and (b) z 2 = 0.5. 



When the DFF is evolved there is a shift towards the 
low Z2 region for z\ fixed at 0.5 (Fig. 4(a)| and towards 
the low zi region for z 2 fixed at 0.5 (Fig. |4(b)[ ). We 
observe that the evolved DFF in Figs. [4] has a larger 
magnitude in the low z\ and z 2 regions, while steadily 
decreasing as z\ and z 2 increase. This is quite different to 
the results shown in Figs. [2j where there is either a large 
contribution for almost all the allowed range of values 
of Z2 (Fig. 2(a) ) or a substantial peak still in the higher 
zi values with the magnitude of the DFF decreasing as 
Zi is decreased. In both those cases, the DFF is largest 
away from the low values of zi and z\. This effect could 
be caused by the down quark, which is produced in both 
fragmentations after the up quark fragments to tt + , being 
an unfavored emission channel for K~ , as opposed to the 
favored emission channel for 7r~ . The favored emission 
channel loses magnitude at higher z\ as the momentum 
scale is increased, while the unfavored emission channels, 
which have no higher z\ peak, increase at lower Z\ due 
to the greater access to the gluon emission channel. 



0.06 
0.05 

(N 
ts) 

J 0.04 

K) 

* 0.03 

+ 

Q 0.02 

Is! 

0.01 



Model Scale 
■ = 0.2 GeV 2 



Evolved Scale 
- =4 GeV 2 



_1_ 



0.2 0.3 




K+K~ 



(c)zi =0.5: z 2 D 



FIG. 5. K + K~ dihadron fragmentation functions for z\ fixed 
to 0.5 at the model scale (Qq — 0.2 GeV 2 , shown by dotted 
red line) and the evolved scale (Q 2 — 4 GeV 2 , shown by solid 
black line) for a fragmenting | (a) | u quark, [(b)] d quark and |(c)| 
s quark. 



C. Q evolution of D„ 



for q = u, d or s 



We now consider D^ +K for q = u (Fig. 5(a) I, 



d (Fig. |5(b)| or s (Fig. |5(c)| ). The q = u and q = sTJFFs 
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both have large peaks in the high z 2 region at the model 
scale since both are favored fragmentation channels in the 



driving function of D 



k+k- 



The first term on the right 



hand side of Eq. |2) produces most of the magnitude of 
the model scale K because the second term equates 
to zero and the intergral term is small. Dff K emerges 
from the second term on the right hand side of Eq. ^ 
because the first term equates to zero and the integral 
term of the DFF is small. The first term on the right 
hand side of Eq. ^ for D*f K contains the elementary 
quark fragmentation function for the fragmentation from 
an up quark to K + as a function of z\, multiplied by 
Df ~0 2 /(l-zi))/(l-zi). For z x fixed to 0.5, this term 
simplifies to a constant multiplied by (z 2 /(\ — z i))- 
However, for DFFs such as K , which emerge from 
the second term on the right hand side of Eq. |2]) , fixing 
Z\ to 0.5 restricts df (z 2 ) to values of z 2 less than 0.5. 
This suppresses the term considerably since the z 2 > 0.5 
region of d^ (z 2 ) is where the function achieves its larger 
values. This is why the u — > K + K~ DFF is larger than 
the s -> K+K- DFF when z x is fixed to 0.5. 

After QCD evolution, the DFFs for fragmenting 
up (Fig. 5(a)) and strange (Fig. |5(c)" ) quarks at z\ = 0.5 
show the shift of the peak value to the lower z 2 region, 
with K having a structure similar to that seen for 



the evolved D£ 77 at z x = 0.5 (Fig. |2(a)| ), while Df K 
has a structure similar to that of the evolved DZ +7! at 



z 2 = 0.5 (Fig. 12(b)). For L>f +K ~ , the model scale plot is 



very small compared to D^ +K and Df +K , since it is 
unfavored for both detected hadrons. When the momen- 



K+K~ 



mcreases m 



turn scale is evolved up to 4 GeV , D: 
the low z 2 region for Z\ fixed to 0.5, because of the effects 
of gluon fragmentation. 



VI. COMPARISON WITH OTHER WORK 

With very little in the way of DFFs from experiments 
being available for comparison, we look to compare our 
results with the work presented in Ref . [75] . We will first 
show that using our code and their parameterized DFFs 
as initial conditions, we do indeed obtain solutions com- 
parable with those presented in Ref. [75J when evolved 
to Q 2 = 109 GeV 2 . We also present our data evolved to 
a range of different scales for -DJ +lr . 

First, we briefly describe the procedures used in 
Ref. [75]. The evolution equations used there are the 



those of Eqs. (15) and (16), with only minor rewriting 



of terms in the equations. For the gluon DFF evolution 
equations, the difference in the equations arises from al- 
ternate definitions of the functions. In Ref. [75], the DFF 
is taken to be identical for the up, down and strange 
quarks, and so the gluon evolution equation term involv- 
ing these functions is written with the function multiplied 
by a factor of 2nt , whereas the DFFs in our approach dif- 
fer and so we sum over each flavor. Similar reasoning is 



used for the other terms in the gluon evolution equation. 
To obtain the initial DFF at Q 2 = 2 GeV 2 , the authors of 
Ref. [75] simulate three million dijet events, distributed 
equally over the number of flavors (rif =3), using JET- 
SET. The resultant DFFs are parameterized by fitting to 
a functional form: 

D( Zl ,z 2 ) =Nz^z^(z 1 +z 2 ) a3 

x (l-z 1 )^(l-z 2 )^(l-z 1 -z 2 )^, (22) 

where N, a±, a 2 , «3, /3i, f3 2 and (3$ are the parame- 
ters fitted by minimizing the logarithm of x 2 . The fit 
describes the JETSET results better at larger values of 
Z\ and z 2 , while not reproducing the results well for 
low values of z\ and z 2 . Values for the parameters are 
provided for the quark and gluon DFFs for momentum 
scales of Q 2 = 2 GeV 2 and Q 2 = 109 GeV 2 . The SFFs 
used are obtained from the parameterization in Ref. [88] . 
The DFFs are QCD evolved from the initial scale of 
Q 2 = 2 GeV 2 , and results are presented for several values 
of Q 2 , including Q 2 = 109 GeV 2 . 

Using the initial parameterized DFFs at Q 2 = 2 GeV 2 , 
in Figs.|6]we present the comparison of the parameterized 
7r+7r~ up quark and gluon DFFs obtained from JETSET 
at Q 2 = 109 GeV 2 (dotted red line) with the evolved 
solutions of Ref. [75] (blue circles). The solutions ob- 
tained using our code on the same initial parameterized 
DFFs (solid black line) and the solution to NJL-jet model 
DFFs evolved to the same momentum scale (solid orange 
line) are shown tocP] We also consider solutions for the 
parameterized DFFs evolved using an altered version of 
our code that treats the QCD evolution of the SFFs with 
the same parameterized evolution as in Ref. |76j (purple 
dot-dashed line), rather than using the evolution equa- 
tions. This serves the purpose of exhibiting how well our 
code reproduces the parameterized JETSET results. 

The results for the NJL-jet model evolved to Q 2 = 
109 GeV 2 are similar to the parameterized JETSET re- 
sults of Ref. |76| for values of z 2 above 0.2 for both the 
up quark (Fig. 6(a)) and gluon (Fig. |6(b)| DFFs. Below 
z 2 = 0.2, our solutions are smaller. Such differences may 
be expected as the parameterization in Ref. |76| overesti- 
mates the JETSET results in the low z\ and z 2 regions, 
and so the NJL-jet model results may well be closer to 
the actual JETSET output. 

We also observe that for the up quark DFF (Fig. |6(a)[ ), 
the solution for the parameterized JETSET input evolved 
using our code produces similar results to the parameter- 
ized solution of the JETSET results at Q 2 = 109 GeV 2 
for values of z 2 greater than approximately 0.1. The 
gluon DFF (Fig. |6(b)[ ) solutions differ only at values of 
z 2 lower than 0.25. In order to understand this differ- 
ence we explored using the parameterized evolution of 



2 These comparisons are at best semi-quantitative as we do not 
know the value of Aqcd used in Ref. 1761 . 
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FIG. 6. 

at Q 2 = 
gluon 



7T 7r~ dihadron fragmentation functions for Zi = 0.5 
= 109 GeV 2 for a fragmenting (a) u quark and 
see text for details. 



(b) 



FIG. 7. 7r + 7r dihadron fragmentation functions 

for the NJL-jet model, for Z\ = 0.5 at Q 2 = 
5 GeV 2 , 20 GeV 2 , 50 GeV 2 and 109 GeV 2 for a frag- 
menting |(a)|u quark and |(b)| gluon. 



the SFFs 88] used by Ref. 76J. This produced an im- 
proved comparison between the parameterized JETSET 
solution and the DFFs obtained through our code. It is 
shown that by employing the parameterized SFF evolu- 
tion we produce results that are similar to the parameter- 
ized JETSET solutions for both the up quark (Fig. 6(a) ) 
and gluon (Fig. |6(b)[ ) for values of z 2 above approxi- 
mately 0.1. 

In Figs.[7J we present the results for the NJL-jet model 
DFFs evolved to a range of Q 2 values: 5 GeV 2 (blue 
dotted line), 20 GeV 2 (black solid line), 50 GeV 2 (green 
dashed line) and 109 GeV 2 (orange solid line). For 
z\ = 0.5, the results show that as Q 2 increases, the DFFs 
appear to gradually reduce. The peak value is also ob- 
served to shift towards lower z 2 values. 



VII. CONCLUSIONS AND OUTLOOKS 

In this article, solutions are presented for dihadron 
fragmentation functions from the NJL-jet model evolved 



to a typical experimental scale of Q 2 = 4 GeV", from 
the model scale of Qq — 0.2 GeV 2 . We first presented 
a brief summary of the integral equations used to ob- 
tain the model scale SFFs and DFFs in Section [TT| Sec- 
tions IIIII and IIVI describe the numerical method used to 
solve the evolution equations for SFFs and DFFs, re- 
spectively. The QCD evolution equations for the SFFs 
and the Fortran code used to solve them was based on 
the method described in Refs. [M l I89H9T] . The method 
used rearranges the evolution equations into non-singlet 
quark and coupled plus-type quark and gluon equations, 
followed by discretizing the variables z and t and con- 
verting the integral terms in to sums over the integration 
variable. The same method is employed to solve the QCD 
evolution equations for the DFFs with the variables z%, 
Z2 and t being discretized. 

Section [V] compares the model scale DFFs with the 
evolved DFFs for 7r + 7r~, ii + K~ and K + K~. In Sec- 
we investigated the evolution of n by com- 



VA 



tion 

paring the model scale and evolved scale DFFs when ei- 
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ther zi (Fig. |2(a)| or z 2 (Fig |2(b)[ ) is eq ual to 0.5. We 
also considered z\ (Fig. 3(a)[ ) or z 2 (Fig. |3(b)[ ) equal to 
0.2. For the fragmentation of the up quark to 7r + 7r~ we 
noted that the up quark was a favored emission channel 
for the 7r + , while the down quark produced after the up 
quark fragments to a tt + is a favored emission channel 
for the 7T~ . The evolved DFF showed a shift in the peak 
value towards the lower z regions, with each plot showing 
the evolved DFF obtaining a larger magnitude than the 
model scale DFF in t he low er z region. 

The focus of Section 
Similar to L>; +7r ~ 
sion channel for 7r 



VB 



is on the evolution of K . 
the up quark is the favored emis- 
+ , however the produced down quark 
is an unfavored emission channel for K~ . The magni- 
tude of the model scale DFF was significantly smaller 
for D? t +K ~ (Figs. [1} than for L>; +7r ~ (Figs. [2} at light- 
cone momentum fractions fixed to 0.5. After evolution, 
D^ +K displayed a similar shift in the peak value to- 
wards the low z region. For z\ — 0.5 (Fig. 4(a)), the 
evolved DFF does not obtain a larger magnitude than 
the model sca le DF F in the lower z-i region, whereas for 
z 2 = 0.5 (Fig. |4(b)| ) the evolved DFF obtains a substan- 
tial increase over the model scale DFF in the lower z\ 
region. This demonstrates the effect evolution has on 
favored and unfavored emission channels. 

Finally, Section |V C| demonstrates the evolution of 
K for q = u, d or s. D*f K has favored fragmen- 
tation channels for both the up quark and strange quark. 
This is observed in the results presented in Figs. [5] where 
both D^ +K and Df +K have large peaks in the upper 
Z2 region. Both D^ +K and Df +K display the shift 
of the peak value to the lower z region that has been 
shown in other favored emission channels at light-cone 
momentum fractions of 0.5. The down quark is an un- 
favored emission for both K + and K~ , and so D^f K 
has a very small magnitude at the model scale. Evolving 
K shows a considerable increase in the lower z 2 re- 
gion, though the magnitude is still much lower than that 



of +K and Df +K at zi = 0.5. 

Evolution of the DFFs has the effect of reducing the 
magnitudes at higher z, resulting in peaks occuring ear- 
lier in the range of z values with a reduced magnitude. If 
the magnitude of the DFF was small at the model scale, 
a significant increase in the magnitude at the low z region 
is observed after evolution. The first of these two effects 
generally occurs for the favored emission channels, where 
the light-cone momentum fraction of the emitted hadron 
is not in the low z region. The second effect, typically 
occurs when the fragmentation channel is unfavored or 
when the emitted hadron carries a small light-cone mo- 
mentum fraction. 

In Section |VI[ , we evolve the parameterized JETSET 
data at 2 Ge\F from Ref. 76 to 109 GeV 2 using our 
code to compare the solutions obtained with the param- 
eterized JETSET data at the same scale for both the up 
quark and gluon fragmenting to 7r + 7r~ (Figs. [6]). We also 
presented solutions for the NJL-jet model up quark and 
gluon DFFs evolved to Q 2 values of 5 GeV 2 , 20 GeV 2 , 
50 GeV 2 and 109 GeV 2 (Figs. [7). The solutions show 
that for z\ — 0.5, the DFFs are reduced as Q 2 increases 
and the peak value shifts towards the lower zi region. 

Extensions of the NJL-jet model for single hadron frag- 
mentation functions such as the inclusion of hadronic 
resonances and their decays [64] and inclusion of the 
transverse momentum dependence [65] have been accom- 
plished using a Monte Carlo framework. These exten- 
sions are possible for DFFs as well, but they are beyond 
the scope of this work are left for the future. 



ACKNOWLEDGMENTS 

This work was supported by the Australian Re- 
search Council through Australian Laureate Fellowship 
FL0992247 (AWT), the ARC Centre of Excellence for 
Particle Physics at the Terascale and by the University 
of Adelaide. 



[1] W. Melnitchouk, J. Speth, and A. W. Thomas, Phys. 

Rev. D 59, 014033 (1998). 
[2] A. W. Thomas, W. Melnitchouk, and F. M. Steffens, 

Phys. Rev. Lett. 85, 2892 (2000). 
[3] A. W. Thomas and W. Weise, The Structure of the Nu- 

cleon (Wiley- VCH, Berlin, Germany, 2001). 
[4] M. Beckmann (HERMES Collaboration), in Testing 

QCD Through Spin Observables in Nuclear Targets, 

edited by D. G Crabb, D. B. Day, and J.-p. Chen (2002), 

pp. 292-298. 

[5] A. Airapetian et al. (HERMES Collaboration), 
Phys.Rev.Lett. 92, 012005 (2004), arXivihep- 
ex/0307064. 

[6] V. Barone, T. Calarco, A. Drago, and M. Simani, 

Phys.Lett. B571, 50 (2003), arXiv:hep-ph/0306225. 
[7] A. Airapetian et al. (HERMES Collaboration), 



Phys.Rev. D71, 012003 (2005), arXiv:hep-ex/0407032. 
[8] A. Airapetian et al. (HERMES Collaboration), 

Phys.Rev. D75, 012007 (2007), arXiv:hep-ex/0609039. 
[9] I. C. Cloet, W. Bentz, and A. W. Thomas, 

Phys.Rev.Lett. 102, 252301 (2009), arXiv: 090 1.3559. 
[10] M. Alekseev et al. (COMPASS Collaboration), 

Phys.Lett. B693, 227 (2010), arXiv:1007.4061. 
[11] F. Gross, G Ramalho, and M. Pena, Phys.Rev. D85, 

093006 (2012), arXiv:1201.6337. 
[12] A. I. Signal and A. W. Thomas, Phys.Lett. B191, 205 

(1987). 

[13] V. Barone, C. Pascaud, and F. Zomer, Eur.Phys.J. C12, 

243 (2000), arXiv:hep-ph/9907512. 
[14] S. Davidson, S. Forte, P. Gambino, N. Rius, and A. Stru- 

mia, JHEP 0202, 037 (2002), arXiv:hep-ph/0112302. 
[15] W. Bentz, I. C. Cloet, J. T. Londergan, and A. W. 



12 



Eur.Phys.J. 
Phys.Rev. D81, 
Phys.Rev. D83, 



114010 



014012 



(2010), 

(2010) , 

(2011) , 



Thomas, Phys.Lett. B693, 462 (2010), arXiv:0908.3198. 
X.-D. Ji, Phys.Rev.Lett. 78, 610 (1997), arXiv:hep- 
ph/9603249. 

X.-D. Ji, Phys.Rev. D55, 7114 (1997), arXiv:hep- 
ph/9609381. 

J. D. Bratt et al. (LHPC Collaboration), Phys.Rev. D82, 
094502 (2010), arXiv:1001.3620. 

G. S. Bali et al. (QCDSF Collaboration), Phys.Rev.Lett. 
108, 222001 (2012), arXiv: 11 12.3354. 

P. Hagler, J.Phys.Conf.Ser. 295, 012009 (2011). 

A. W. Thomas, A. Casey, and H. H. Matevosyan, 

Int.J.Mod.Phys. A25, 4149 (2010). 

A. W. Thomas, Phys.Rev.Lett. 101, 102003 (2008), 
arXiv:0803.2775. 

F. Myhrer and A. W. Thomas, Phys.Lett. B663, 302 
(2008), arXiv:0709.4067. 

S. D. Bass and A. W. Thomas, Phys.Lett. B684, 216 

(2010) , arXiv:0912.1765. 

S. D. Bass, A. Casey, and A. W. Thomas, Phys.Rev. 
C83, 038202 (2011), arXiv:1110.5160. 
M. Wakamatsu, Eur.Phys.J. A44, 297 
arXiv:0908.0972. 
M. Wakamatsu, 
arXiv:1004.0268. 
M. Wakamatsu, 
arXiv:1007.5355. 

A. Adare et al. (PHENIX Collaboration), Phys.Rev.Lett. 

103, 012003 (2009), arXiv:0810.0694. 

A. Bacchetta, D. Boer, M. Diehl, and P. J. Mulders, 

JHEP 0808, 023 (2008), arXiv:0803.0227. 

V. Barone, A. Drago, and P. G. Ratcliffe, Phys.Rept. 

359, 1 (2002), arXiv:hep-ph/0104283. 

E. Leader, A. V. Sidorov, and D. B. Stamenov (2010), 

arXiv:1012.5033. 

E. Leader, A. V. Sidorov, and D. B. Stamenov, Phys.Rev. 

D82, 114018 (2010), arXiv:1010.0574. 

E. Leader, Phys.Rev. D83, 096012 (2011), 

arXiv:1101.5956. 

E. Leader, A. V. Sidorov, and D. B. Stamenov, 
J.Phys.Conf.Ser. 295, 012054 (2011). 

W. Vogelsang, Nucl.Phys. A827, HOC (2009). 

F. A. Ceccopieri, Czech. J.Phys. 56, F175 (2006), 
arXiv:hep-ph/0610074. 

A. Bacchetta, M. Diehl, K. Goeke, A. Metz, P. J. 
Mulders, et al., JHEP 0702, 093 (2007), arXiv:hep- 
ph/0611265. 

I. C. Cloet, W. Bentz, and A. W. Thomas, Phys.Lett. 
B659, 214 (2008), arXiv:0708.3246. 

L. P. Gamberg, A. Mukherjee, and P. J. Mulders, 

Phys.Rev. D83, 071503 (2011), arXiv: 1010.4556. 

L. Gamberg and M. Schlegel, AIP Conf.Proc. 1374, 309 

(2011) , arXiv:1012.3395. 

H. Wollny (COMPASS Collaboration), in Exclusive Re- 
actions at High Momentum Transfer IV, edited by 
A. Radyushkin (World Scientific, Singapore, 2011), pp. 
303-311. 

M. Ansclmino, H. Avakian, D. Boer, F. Bradamante, 
M. Burkardt, et al., Eur.Phys.J. A47, 35 (2011), 
arXiv:1101.4199. 

M. Anselmino, M. Boglione, U. D'Alesio, S. Melis, 
F. Murgia, et al., Phys.Rev. D83, 114019 (2011), 
arXiv:1101.1011. 

S. M. Aybat, A. Prokudin, and T. C. Rogers (2011), 
arXiv:1112.4423. 



[46 
[47 

[48; 

[49 

[50; 

[51 
[52 

[53 

[54 

[55 

[56 

[57; 

[58 
[59 
[60; 
[61 
[62 

[63 
[64 
[65 

[66 
[67 
[68; 
[69 
[70 
[71 
[72 

[73 

[74 



L. Gamberg, D. Boer, B. Musch, and A. Prokudin, AIP 
Conf.Proc. 1418, 72 (2011), arXiv:1111.0603. 
A. Bacchetta and M. Radici, Phys.Rev.Lett. 107, 212001 
(2011), arXiv:1107.5755. 

K. Goeke, M. V. Polyakov, and M. Vanderhaeghen, 
Prog.Part. Nucl.Phys. 47, 401 (2001), arXivihep- 
ph/0106012. 

M. Burkardt, Int.J.Mod.Phys. A18, 173 (2003), 
arXiv:hep-ph /0207047. 

M. Diehl, Phys.Rept. 388, 41 (2003), arXiv:hep- 
ph/0307382. 

X. Ji, Ann.Rev.Nucl.Part.Sci. 54, 413 (2004). 

A. V. Belitsky and A. V. Radyushkin, Phys.Rept. 418, 

1 (2005), arXiv:hep-ph/0504030. 

S. Boffi and B. Pasquini, Riv.Nuovo Cim. 30, 387 (2007), 
arXiv:071 1.2625. 

M. Radici, II Nuovo Cimento C 035, 02 69 (2012), 
arXiv:1111.3383. 

M. Hirai, S. Kumano, T.-H. Nagai, and K. Sudoh, 
Phys.Rev. D75, 094009 (2007), arXiv:hep-ph/0702250. 
D. de Florian, R. Sassot, and M. Stratmann, Phys.Rev. 
D75, 114010 (2007), arXiv:hep-ph/0703242. 
O. Biebel, P. Nason, and B. R. Webber, Tech. Rep. 
hep-ph/0109282. BICOCCA-FT-2001-20. CAVENDISH- 
HEP-2001-12, Milano-Bicocca Univ. Phys. Dept., Milano 
(2001). 

D. de Florian and L. Vanni, Phys.Lett. B578, 139 (2004), 
arXiv:hep-ph/0310196. 

M. Hirai, S. Kumano, T.-H. Nagai, and K. Sudoh, AIP 
Conf.Proc. 915, 749 (2007), arXiv:hep-ph/0612009. 
M. Hirai, S. Kumano, T.-H. Nagai, M. Oka, and K. Sudoh 
(2007), arXiv:0709.2457. 

P. Francisconi, Nucl.Phys.Proc.Suppl. 207-208, 129 
(2010). 

T. Ito, W. Bentz, I. C. Cloet, A. W. Thomas, 
and K. Yazaki, Phys.Rev. D80, 074008 (2009), 
arXiv:0906.5362. 

H. H. Matevosyan, A. W. Thomas, and W. Bentz, 

Phys.Rev. D83, 074003 (2011), arXiv:1011.1052. 

H. H. Matevosyan, A. W. Thomas, and W. Bentz, 

Phys.Rev. D83, 114010 (2011), arXiv:1103.3085. 

H. H. Matevosyan, W. Bentz, I. C. Cloet, and 

A. W. Thomas, Phys.Rev. D85, 014021 (2012), 

arXiv:1111.1740. 

R. D. Field and R. P. Feynman, Nucl.Phys. B136, 1 
(1978). 

H. H. Matevosyan, A. W. Thomas, and W. Bentz (2012), 
arXiv: 1205.5813. 

H. H. Matevosyan, A. W. Thomas, and W. Bentz (2012), 
arXiv: 1207.0812. 

H. H. Matevosyan, A. W. Thomas, and W. Bentz (2012), 
arXiv: 1207. 1433. 

A. Bacchetta and M. Radici, Phys.Rev. D74, 114007 
(2006), arXiv:hep-ph/0608037. 

J. Zhou and A. Metz, Phys.Rev.Lett. 106, 172001 (2011), 
arXiv:1101.3273. 

T. Sjostrand, P. Eden, C. Friberg, L. Lonnblad, 
G. Miu, et al., Comput.Phys.Commun. 135, 238 (2001), 
arXiv:hep-ph/0010017. 

P. Liebing, Ph.D. thesis, Universitat Hamburg (2004), 
[DESY Report No. DESY-THESIS-2004-036, 2004 (un- 
published)]; http://cdsweb.cern.ch/record/808379. 
A. Casey, H. H. Matevosyan, and A. W. Thomas, Phys. 
Rev. D 85, 114049 (2012). 



13 



[75] A. Majumder and X.-N. Wang, Phys.Rev. D70, 014007 

(2004) , arXiv:hep-ph/0402245. 

[76] A. Majumder and X.-N. Wang, Phys.Rev. D72, 034007 

(2005) , arXiv:hep-ph/0411174. 

[77] L. Grigoryan (2008), arXiv:0809.0281. 

[78] S. Boffi, R. Jakob, M. Radici, and A. Bianconi, in Per- 
spectives in Hadronic Physics: Proceedings, 2nd Inter- 
national Conference, Trieste, Italy, May 10 -14, edited 
by S. Boffi, C. Ciofi degli Atti, and M. Giannini (World 
Scientific, Singapore, 2000), pp. 467-478, arXiv:hep- 
ph/9907374. 

[79] M. Radici, in Physics With A High Luminosity Polarized 
Electron Ion Collider: Proceedings., Bloomington, Indi- 
ana, USA, 8-11 April 1999, edited by L. C. Bland, J. T. 
Londergan, and A. P. Szczepaniak (World Scientific, Sin- 
gapore, 2000), pp. 113-122, arXiv:hep-ph/9906217. 

[80] M. Radici, in Gerasimov-Drell-Hearn Sum Rule and the 
Spin Structure of the Nucleon: Proceedings, 2nd Inter- 
national Symposium, GDH 2002, Genova, Italy, July 3 
- 6, 2002, edited by M. A. Anghinolfi, M. Battaglieri, 
and R. de Vita (World Scientific, Hackensack, NJ, 2003, 
2002), pp. 397-402, arXiv:hep-ph/0209348. 

[81] A. Bacchetta and M. Radici, Phys.Rev. D70, 094032 
(2004), arXiv:hep-ph/0409174. 

[82] J. She, Y. Huang, V. Barone, and B.-Q. Ma, Phys.Rev. 



D77, 014035 (2008), arXiv:0711.0817. 

[83] A. Bacchetta, A. Courtoy, and M. Radici, Phys.Rev. Lett. 
107, 012001 (2011), arXiv:1104.3855. 

[84] A. Courtoy, A. Bacchetta, and M. Radici, in proceed- 
ings of XIX International Workshop on Deep-Inelastic 
Scattering and Related Subjects (DIS 2011), April 11-15, 
2011, Newport News, VA USA (2011), arXiv:1106.5897. 

[85] A. Courtoy, A. Bacchetta, M. Radici, and A. Bianconi, 
Phys. Rev. D 85, 114023 (2012). 

[86] M. Hirai and S. Kumano, Comput.Phys.Commun. 183, 
1002 (2012), arXiv:1106.1553. 

[87] F. A. Ceccopieri, M. Radici, and A. Bacchetta, 
Phys.Lett. B650, 81 (2007), arXiv:hep-ph/0703265. 

[88] J. Binnewies, B. A. Kniehl, and G. Kramer, Phys. Rev. 
D 52, 4947 (1995). 

[89] M. Miyama and S. Kumano, Comput.Phys.Commun. 94, 
185 (1996), arXiv:hep-ph/9508246. 

[90] M. Hirai, S. Kumano, and M. Miyama, in High-energy 
spin physics. Proceedings, 12th International Symposium, 
SPIN 96, Amsterdam, Netherlands, September 10-14, 
1996., edited by C. W. de Jager, T. J. Ketel, P. J. Mul- 
ders, J. E. J. Oberski, and M. Oskam-Tamboezer (1996), 
pp. 410-412, arXiv:hep-ph/9610521. 

[91] M. Hirai, S. Kumano, and M. Miyama, Com- 
put.Phys.Commun. 108, 38 (1998), arXiv:hep- 
ph/9707220. 



